Competition between two- and three-sublattice ordering for 5 = 1 spins on the square lattice 
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We provide strong evidence that the S = 1 bilinear-biquadratic Heisenberg model with nearest-neighbor 
interactions on the square lattice possesses an extended three-sublattice phase induced by quantum fluctuations 
for sufficiently large biquadratic interactions, in spite of the bipartite nature of the lattice. The argumentation 
relies on exact diagonalizations of finite clusters and on a semiclassical treatment of quantum fluctuations within 
linear flavor-wave theory. In zero field, this three-sublattice phase is purely quadrupolar, and upon increasing 
the field it replaces most of the plateau at 1/2 that is predicted by the classical theory. 
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On bipartite lattices, classical magnetic models with 
nearest-neighbor interactions usually display very simple 
types of order. Indeed, once the configuration that minimizes 
the energy of a pair of sites is known, the total energy can be 
minimized by extending the solution to the two sublattices of 
the bipartite lattice. If both sites of the pair are in the same 
state, the resulting configuration is uniform; otherwise it cor- 
responds to a two-sublattice state. This is true for Ising, XY, 
and Heisenberg models, with a uniform ferromagnetic ground 
state for negative coupling and a two-sublattice Neel ground 
state for positive coupling. This reasoning can be extended to 
quantum models by considering product wavefunctions (see 
below), and if quantum fluctuations around this "classical" 
solution do not destroy the order, the ground state can be ex- 
pected to be uniform or Neel-Iike. This is indeed known to 
be the case for the Heisenberg model for any spin. One might 
naively assume that this simple argument also applies to the 
generic spin-1 SU(2)-invariant model, the bilinear-biquadratic 
Hamiltonian 



(1) 



since although the interaction includes a biquadratic term, it 
is still limited to nearest neighbors. 

However, this picture is challenged by two observations. 
First of all, it is well established by now that in the regime 
< Ji < J2, the correlations of the ID case are algebraic 
with period 3.'"^ This indicates that weakly coupled chains 
will be driven to a long-range ordered state that respects this 
periodicity, and this might also occur in the isotropic 2D case. 
Second, our recent study of the SU(3) Heisenberg model on 
the square lattice'* has shown that three-flavor stripe order is 
stabilized, and this model is equivalent to the model of Eq. (1) 
with Ji = J2 in zero field. This raises a natural question: 
Is the three-sublattice order a consequence of the enhanced 
SU(3) symmetry, or does the model of Eq. (1) actually possess 
a phase with this symmetry? 

In this Rapid Communication, we show that quantum fluc- 
tuations actually stabilize a three-sublattice phase in an un- 
expectedly large portion of the phase diagram, including a 



parameter range in zero field where chains are gapped with 
incommensurate short-range correlations,^ and a finite-field 
region where classical spins form a two-sublattice magneti- 
zation plateau at 1/2. 

In order to set the stage, let us begin with a review of the 
classical (Hartree) solution, a variational approach based on a 
product wave function of the form 



N 



(2) 



where N is the number of lattice sites. The minimization of 
the energy with respect to the local wave functions \ipi) was 
achieved by Papanicolaou some time ago.^ Given the impor- 
tance of the results for the discussion of quantum effects, let 
us present the solution in a slightly more compact way than 
Ref. 6 does. In order to discuss general spin-1 states, it is 
convenient to introduce the time-reversal invariant basis 
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\z) = -i\0). (3) 



A general state |d) = da\a) can be described by a com- 
plex vector d = u + tv^. Without loss of generality, one can 
choose u and v in such a way that u • v = 0, and the normal- 
ization imposes u'^ + ~ 1. If u = v, the state is purely 
magnetic and (S) = 2 u x v. If u = or u = 0, the state 
is purely quadrupolar with a director along the nonzero com- 
ponent u or v: It is obtained by an SU(2) rotation from |0) 
and it displays spin fluctuations in a plane perpendicular to its 
director. Following the general strategy for bipartite lattices, 
we first minimize the variational energy of a pair of sites: It is 
given by 



E = Ji|d, • d*|2 + (J2 - Ji)|d, • d,f + J2 



(4) 



d*|2 and 



This energy is minimum for configurations with |d 
\di ■ dj p equal to or 1 (it turns out that the extreme values 
can be attained simultaneously) depending on whether Ji and 
J2 — Ji are positive or negative. This leads to the following 
solutions: 
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FIG. 1: (Color online) Schematic phase diagram of the S = I 
bilinear-biquadratic model on the square lattice. The inner circle 
shows the variational, the outer the numerical results. FM and AFM 
are ferro- and antiferromagnetic, FQ is ferroquadrupolar, while SO 
stands for semiordered. In the outer circle, we have shown the three- 
sublattice ordered antiferromagnetic and antiferroquadrupolar phases 
between -d ~ 0.2n and tt/2. 
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FIG. 2: (Color online) (a) The variational phase diagram in a mag- 
netic field. AFM2 denotes a two-sublattice ordered Neel phase, while 
FM stands for ferromagnet. Note the presence of a magnetization 
plateau at 1/2 which is separated from AFM2 by a supersolid phase, 
(b) Schematic plot of the anticipated phase diagram based on ex- 
act diagonalization calculations. AFM3 (AFQ3) is a three-sublattice 
ordered antiferromagnetic (antiferroquadrupolar) phase. The upper 
boundary of the 1/2 plateau is a result of finite-size scaling. The 
dashed line is estimated within flavor- wave theory (see text). 



i) Ji > 0, J2 — Ji < 0; dj = d* to maximize |di-djp, and 
Ui = Vi to satisfy • d* = uf — vf ~ 0. The states are purely 
magnetic and the spins are antiparallel. Note that this case 
includes the simple antiferromagnetic Heisenberg coupling. 

ii) Ji < 0, J2 — Ji < 0: dj = d* to maximize \di ■ dj p, 
and = or = to maximize |dj • d*p = {uf — vf)^. 
The states are purely quadrupolar with parallel directors. 

iii) Ji < 0, J2 — Ji > 0: dj = di to maximize |di • djp, 
and Ui ~ Vi to satisfy d^ • dj = — vf = 0. The states are 
purely magnetic with parallel spins. This case includes the 
simple ferromagnetic coupling. 

iv) Ji > 0, J2 — Ji > 0: dj • dj = d; • d* = implies 
Ui • Uj — Ui ■ Vj — Vi ■ Uj — Vi ■ Vj = 0. Since • = 
Uj ■ Vj =0 by definition, the four vectors u^, v^, Uj, and Vj 
must be pairwise orthogonal. This implies that at least one 
of them vanishes; i.e., one state is a pure quadrupole. The 
other state is only constrained by the condition that its u and 
V vectors must be perpendicular to the director of the first 
state. It can be purely magnetic with a moment along this 
director, purely quadrupolar with a director perpendicular to 
this director, or of mixed character. Unlike the other cases, 
where the ground state is unique up to a global SU(2) rotation, 
the ground state is highly degenerate in this case. 

The above analysis leads to the zero-field classical phase 
diagram shown in Fig. 1 (inner circle), with the standard no- 
tation Ji = Jcos?9 and J2 = Jsind. It consists of four 
phases, clockwise from the Heisenberg point: antiferromag- 
netic (AFM), ferroquadrupolar (FQ), ferromagnetic (FM), and 
semiordered (SO) in the terminology of Papanicolaou.^ The 
last phase is highly degenerate. For instance, if all sites of 
one sublattice are in the same quadrupolar state, the sites 
of the other sublattice can choose their state (purely mag- 
netic/quadrupolar or mixed) independently of each other An- 
other family consists of all coverings of the square lattice us- 
ing three quadrupoles with mutually perpendicular directors. 



with a residual entropy equal to that of the 3-state Potts model 
on the square lattice. 

In a magnetic field, this degeneracy is lifted: The two- 
sublattice configuration with identical quadrupoles on one 
sublattice and purely magnetic states with a moment along the 
field on the other sublattice is the unique ground state. This 
gives rise to a magnetization plateau at 1/2 (see Fig. 2). Upon 
leaving the plateau, the system adopts a two-sublattice con- 
figuration with mixed states: At first, these states feature dif- 
ferent polarizations along the field (a kind of spin supersolid), 
and subsequently a canted antiferromagnet is stabilized. 

In summary, up to the degeneracy of the ground state in 
zero field for 7r/4 < < tt/2, the classical phase diagram 
consists of uniform or two-sublattice phases, with an abrupt 
transition into a 1/2 plateau for an infinitesimal field in this 
parameter range. As we shall now demonstrate, this picture is 
considerably affected by quantum fluctuations. 

Numerical approach. For J2 < 0, there is no sign prob- 
lem for quantum Monte Carlo simulations, and Harada and 
Kawashima** have shown that the sequence of phases (AFM, 
FQ and FM) as well as the location of the transitions predicted 
by the classical approach are preserved by quantum fluctua- 
tions. In the following, we concentrate on the parameter range 
Ji, J2 > (0 < I? < 7r/2). In order to get an insight into the 
nature of the zero-field phase, we have calculated the structure 
factors cxp[ik • Rj](C(0) • C(Rj)), where C(Rj) is the 
spin or quadrupole operator^ at site Rj, by exact diagonaliza- 
tion for different wave vectors k as a function of ^ (Fig. 3). 
For -d ~ (Heisenberg model), the structure factor is the 
largest (and grows with the system size) at (tt, tt), as we ex- 
pect for a two-sublattice ordered Neel antiferromagnet. As we 
turn on the biquadratic coupling, the structure factor gradually 
decreases, and for i? > O.IQtt, i.e., long before the classical 
transition point at 7r/4 is reached, three-sublattice correlations 
take over, with a structure factor that peaks at {2tt/3, 2tt/3). 
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Spin-spin correlations prevail for i? < 7r/4, while quadrupolar 
correlations become dominant in the region where J2 > Ji. A 
three-sublattice ordering is also suggested by the peculiar de- 
pendence of the ground-state energy on the number of sites: 
We consistently get lower energies for clusters that are mul- 
tiples of three. Furthermore, for > 7r/4, the low-energy 
spectrum of the 18-site cluster reveals an Anderson tower"''^ 
(Fig. 4) that corresponds to two copies of the Anderson tower 
of the three-sublattice antiferroquadrupolar state on the trian- 
gular lattice. '^"'^ The two copies refer to the Z2 degeneracy 
of the state, i.e., to the orientation of the stripes. We expect 
a finite-temperature ordering transition of Ising type to select 
the orientation of the stripes, while the antiferroquadrupolar 
order is only established at T = 0. 

We complete our numerical analysis by investigating the 
magnetization process in the range 7r/8 < -Q < 7r/2 (Fig. 5). 
In contrast to the classical solution, the magnetization does not 
jump directly to 1/2 for 7r/4 < ?? < 7r/2: There is an interme- 
diate phase at low field in which the system polarizes progres- 
sively and retains the three-sublattice order of the zero-field 
ground state. In order to see whether the plateau at 1/2 re- 
mains in the thermodynamic limit, we have attempted a finite- 
size scaling of its critical fields for 16, 18, and 20 sites. For 
the upper critical field, a linear finite-size scaling is possible 
for -d not too close to 7r/4, resulting in the upper boundary of 
Fig. 2(b). For the lower critical field, however, no meaning- 
ful finite-size scaling could be performed, presumably due to 
the incompatibility of the 16- and 20-site clusters with three- 
sublattice long-range order 

Semiclassical approach. The numerical results are con- 
firmed and complemented by a semiclassical analysis that 
includes zero-point fluctuations at the level of linear flavor- 
wave theory.''-''"''^ In zero field for tt/4 < < tt/2, we 
have compared the energy of all two-sublattice ordered classi- 
cal ground states and the three-sublattice antiferroquadrupo- 
lar state, and the latter state is indeed stabilized for all 1!). 
For d = IT /A, this state is considerably lower in energy than 
the two-sublattice antiferromagnetic state, and even though 
we cannot follow the three-sublattice state below d = tt/A, 
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FIG. 3: (Color online) Spin and quadrupole structure factor from ex- 
act diagonalizations of square clusters with periodic boundary condi- 
tions. The geometry of the clusters can be found in Fig. 72 of Ref. 10. 
Open (closed) symbols stand for the spin (quadrupole) structure fac- 
tor. 
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FIG. 4: (Color online) Anderson tower for the 18-site cluster show- 
ing the two-sublattice (a) and the three-sublattice order (b). The 
dashed line is a guide to the eye. The energy levels are character- 
ized by their momenta and irreducible representations. 
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FIG. 5: (Color online) Magnetization curves of finite clusters for 
i9 = Stt/S. Dashed lines represent variational magnetization curves 
calculated for two-sublattice and three-sublattice ordered states. The 
3D plot shows the magnetization process of the 18-site cluster as a 
function of i9. 



since it is no longer a classical ground state and imaginary 
frequencies appear in its flavor-wave spectrum, we may argue 
by continuity that three-sublattice order will be favored over 
two-sublattice antiferromagnetic order in a finite window be- 
low ?? = 7r/4, in agreement with exact diagonalizations. In 
a finite field for tt/A < -d < tt/2, the 1/2-plateau state is the 
unique classical ground state and again we cannot calculate 
the energy of the three-sublattice state. However, for small 
fields, the correction to the energy is dominated by the classi- 
cal energy gain of the 1/2-plateau state, which is linear in field 
and equal to h/2, since all other corrections are quadratic. 
Therefore, the boundary between the three-sublattice phase 
and the 1/2 plateau is given hy h = 2(e2 — £3), where £3 
(resp. £2) is the zero-point correction to the energy of the 



4 



(a) 


O-0.26II 




h^J=0.54 - 






\ 


/ 




(ii/2,n/2) (0,0) 



(ii,0)(jt/2,ii/2) (0,0) 

k 



(ii,0) (ii/2,ii/2) (0,0) 



(11,0) 



FIG. 6: (Color online) Dynamical structure factors for the plateau 
state in the reduced Brillouin zone for three different values of -d at 
the upper critical magnetic field where the gap closes. The color 
code goes from yellow for spin waves to red for quadrupolar waves 
according to the dominant character of the mode. The dispersing 
bands have finite 5"" = 5" and Q"'' = (J^"" matrix elements. The 
flat dispersion at (xj = 2/i is a pure quadrupolar wave with Q^^ = 

2 _ 2 

^ fluctuations. The spectra are displayed at the critical field to 
show that it is a partially magnetic mode that softens at the transition. 



three-sublattice antiferroquadrapolar (resp. 1/2 -plateau) state 
calculated in zero field. This boundary is shown in Fig. 2(b). 
In the I? — )- 7r/2 limit, the critical field goes to zero and the 
approximation is quantitatively reliable. Now, to leading or- 
der in 5 = 7r/2 — t?, the zero-point corrections are given by 
£3 = -8(5^/^/37r+0((52) and = In (5; therefore the criti- 
cal field scales as /i cx 8'^!'^. This implies that the lower bound- 
ary of the plateau has a vanishing slope in the d ^ n jl limit. 
Since the finite-size scaling of the upper critical field deduced 
from exact diagonalizations is very good in the same limit and 
yields a significant slope, one can safely conclude that a 1/2 
plateau is indeed present in the phase diagram. However, the 
fate of this plateau upon approaching ?9 = 7r/4 cannot be de- 
duced from this argument, and the corresponding region of the 
phase diagram of Fig. 2 is just a plausible scenario. Note that 
both the semiclassical calculation and exact diagonalizations 
are consistent with a magnetization jump at the lower critical 
field of the 1/2 plateau. 

Finally, let us comment on the nature of the elementary ex- 
citations. It is best revealed by the dynamical structure factors 

C(k,c.) = ^|(HQ(k)|G)p5(c.-^,), (5) 



where Cj (k) can be the Fourier transform of a spin or a 
quadrupole operator on sublattice j, the ground state is de- 
noted by |G), and \v) are excited states with energy w^. In 
general, the excitations are of mixed character; i.e., they ap- 
pear with measurable weight both in the spin and quadrupole 
dynamical structure factors, with on average a stronger spin 
resp. quadrupolar character depending on the type of order, 
on the band, and on the wave vector For the 1/2-plateau 
phase, however, the flavor-wave analysis predicts that among 
the four bands, one is dispersionless and of pure quadrupolar 
character for all as shown in Fig. 6. This mode would be 
completely absent in neutron scattering, but could be detected, 
e.g., in Raman scattering.'^ 

The interpretation of the spectrum is particularly transpar- 
ent for = 7r/2, where the variational plateau state becomes 
an exact ground state. Restricting ourselves to single-particle 
states, we find that the \z) |1) (with AS"^ = +1) and 
|1) |1) (AS"" = —2) excitations are localized, forming flat 
bands. The |z) ^ |I) and |1) ^ |z) excitations (A5^ = -1) 
on the two sublattices hybridize and form two bands with dis- 
persion [4 ± 2(cosfci, + cosfcj,)] J. These excitations are re- 
covered in the flavor-wave calculations in the d ^ 7^/2 limit, 
and can be recognized in Fig. 6(c). Away from d ~ 7r/2, 
the \z) —J- |1) excitation acquires a dispersion, while the 
|1) |i) quadrupolar excitation remains dispersionless. 

Conclusions. In summary, the evidence presented here 
strongly suggests that the bilinear-biquadratic Heisenberg 
model on the square lattice features a phase with three- 
sublattice order, although the lattice is bipartite and the cou- 
pling is limited to nearest neighbors. The extent of this phase 
is quite surprising in view of the classical phase diagram; It 
replaces a large portion of the 1/2 plateau in a field, and it 
extends to the region where the exchange is predominantly bi- 
linear The 1/2 plateau is also quite interesting, with a purely 
quadrupolar dispersionless excitation. It is our hope that the 
present work will further motivate the search for S = 1 quan- 
tum magnets with significant positive biquadratic interactions. 
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